# Generate Figure 4
set.seed(102112516) 

N <- 15000
m <- N

# Generate the appropriate lognormals according to the paper.
X1 <- mvrnorm(N,c(0,0),matrix(c(1.2,0,0,1.2),nrow=2))
Y1 <- exp(X1)

X2 <- mvrnorm(N,c(0,0),matrix(c(1,0.9,0.9,1),nrow=2))
Y2 <- exp(X2)

ilf1 <- ILF(Y1,rep(1/N,N),vquantile(Y1),m)
ilf2 <- ILF(Y2,rep(1/N,N),vquantile(Y2),m)

# Calculate the coordinates of the corners of the identical distribution
# There are more alphas than necessary.
alpha <- c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.75,0.8,0.9,0.95,0.99)
egal <- rep(0,12)
for (k in 1:12){
  # This formula is given in paper section 1.2.2 example 1, set d = 2
  egal[k] <- uniroot(function(x){x*(1-log(x)) - alpha[k]},c(0.001,1))$root
}

# Generate grid
sequence <- rep(seq(from = 0, to = 1, by = 1/(nrow(ilf1$estimate)-1)))
s <- 0
t <- 0
x <- rep(0,length(sequence)^2)
y <- x
for (i in 1:length(sequence)){
  s <- s+1
  k <- length(sequence)*(s - 1) + 1
  j <- length(sequence)*s
  x[k:j] <- rep(sequence[s],length(sequence))
  y[k:j] <- sequence
  
}

lr <- data.frame(x,y,c(ilf1$estimate),c(ilf2$estimate))
names(lr) <- c("x","y","X_1","X_2")
lr <- reshape2::melt(lr, id.var = c("x","y"))

levels <- c(0.2,0.5,0.75,0.9,0.95)
sz <- 1

lcontours <- ggplot(data=lr,aes(x=x,y=y,z=value)) + 
  geom_contour(aes(colour = variable),breaks = levels,size =sz) +
  annotate("label", x=0.9, y=0.41, label= "0.95", family = 'serif',color ="black",fill = "white", size = 5.5) +
  annotate("label", x=0.9, y=0.28, label= "0.90", family = 'serif',color ="black",fill = "white", size = 5.5) +
  annotate("label", x=0.9, y=0.13, label= "0.75", family = 'serif',color ="black",fill = "white", size = 5.5) +
  annotate("label", x=0.9, y=0.05, label= "0.50", family = 'serif',color ="black",fill = "white", size = 5.5) +
  annotate("label", x=0.8, y=0.0, label= "0.20", family = 'serif',color ="black",fill = "white", size = 5.5) +
  annotate("segment",x = 0, y = 0, xend = 0.85, yend = 0.85,linetype = 3, size=sz)+
  scale_color_manual(labels = c(expression(paste(sigma^2," = ", "1.2",sep = "")),expression(paste(sigma^2," = ", "1.0",sep = ""))), values=c("X_1" = "black","X_2"="darkgray")) +
  coord_cartesian(xlim = c(0, 1),ylim = c(0, 1)) +
  theme(plot.title = element_text(hjust = 0.5)) + xlab("")+ ylab("")+ 
  guides(col= guide_legend(title= "Scale"))+
  theme(panel.border = element_rect(colour = "black", fill=NA, size=0.6),
        panel.background = element_rect(fill = 'transparent'),
        panel.grid = element_blank(),
        legend.position = c(0.85,0.85),
        legend.key = element_rect(fill = 'transparent'),
        text = element_text(family = "serif",size=22))+
  scale_y_continuous(breaks = seq(0,1,0.2),position = "right",labels = seq(0,1,0.2)) +
  scale_x_continuous(breaks = seq(0,1,0.2),position = "top",labels = seq(0,1,0.2)) 
ggsave(plot = lcontours, width = 6, height = 6, filename = "fig4.png")


